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Abstract 

Cooperative Greedy Pursuit Strategies are considered for approximating a signal partition 
subjected to a global constraint on sparsity. The approach aims at producing a high quality 
sparse approximation of the whole signal, using highly coherent redundant dictionaries. The co¬ 
operation takes place by ranking the partition units for their sequential stepwise approximation, 
and is realized by means of i)forward steps for the upgrading of an approximation and/or ii) 
backward steps for the corresponding downgrading. The advantage of the strategy is illustrated 
by approximation of music signals using redundant trigonometric dictionaries. In addition to 
rendering stunning improvements in sparsity with respect to the concomitant trigonometric 
basis, these dictionaries enable a fast implementation of the approach via the Fast Fourier 
Transform. 


Keywords: Hierarchized Block Wise Pursuit Strategies, Optimized Orthogonal Matching Pursuit, 
Sparse Representation of Music Signals. 

1 Introduction 

Compressible signals, such as audio and vision data, are characterized by samples containing a good 
deal of redundancy. Consequently, if properly processed, the signal information content can be 
accessed from a reduced data set. Transformations for data reduction are said to produce a sparse 
representation of a signal if they can accurately reproduce the information the signal conveys, with 
significantly less points that those by which the original signal is given. Most popular transformations 
for signal processing (e.g. Discrete Cosine Transform and Discrete Wavelet Transform) do not modify 
the signal size. A sparse representation is achieved, a posteriori, by disregarding the least relevant 
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points in the transformed domain. Nonetheless, much higher sparsity in a signal representation 
may be achieved by allowing for the expansion of the transformed domain, thereby leaving room 
for dedicated transformations adapted to the particular signal. Such a framework involves a large 
redundant set called ‘dictionary’. The aim is to represent a signal as a superposition of, say K, 
dictionary’s elements, which are called ‘atoms’, with K much less than the signal size. 

Given a redundant dictionary, the problem of hnding the sparsest approximation of a signal, up 
to some predetermined error, is an NP-hard problem [1]. Consequently, the interest in practical 
applications lies in the hnding of those solutions which, being sparse in relation to other possible 
representations, are also easy and fast to construct. Sparse practical solutions can be found by 
what are known as Pursuit Strategies. Here the discrimination between two broad categories is 
in order: i)The Basis Pursuit based approaches, which endeavor to obtain a sparse solution by 
minimization of the 1-norm [2]. ii) Greedy algorithms which look for a sparse solution by stepwise 
selection of dictionary’s atoms. Practical greedy algorithms, which originated as regression techniques 
in statistics [3], have been popularized in signal processing applications as Matching Pursuit (MP) [4] 
and Orthogonal Matching Pursuit (OMP) [5] methods. The approach, which in principle consider 
the stepwise selection of single atoms, has been extended to multiple atom selection [6]. Dedicated 
algorithms such as Stagewise Orthogonal Matching Pursuit [7], Gompressive Sampling Matching 
Pursuit [8], and Regularized Orthogonal Matching Pursuit [9] are known to be effective within the 
context of the emerging theory of sampling, called compressive sensing/sampling. This theory asserts 
that sparsity of a representation may also lead to more economical data collection [10-14]. In that 
context the reconstruction problem is of a very particular nature, though: It is assumed that a signal 
is sparse in an orthogonal basis and the goal is to reconstruct the signal from a reduced number 
of measures. On the contrary, in this Gommunication we address the traditional representation 
matter: the signal is assumed to be completely given by its samples. The aim is to produce a high 
quality approximation of all those samples, as a iP-term superposition of atoms belonging to a highly 
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coherent dictionary. In this case, minimization of the 1-norm is not effective and step wise greedy 
selection of single atoms benefits sparsity results. 

In practice, when trying to approximate real life signals using a redundant dictionary, there is 
a need to approximate by partitioning. While this requirement normally comes from storage and 
computational complexity demands, it does not represent a disadvantage. On the contrary, to capture 
local properties of signals such as images, or music, non-local orthogonal transforms are also applied 
on partitions to obtain superior results. The broadly used image compression standard JPEG and 
the music compression standard MP3, for instance, operate by partitioning as a hrst step in the 
compression process. 

The central aim of this paper is to tailor pursuit algorithms for approximation by partitioning 
without signihcant increment in computational complexity, in comparison with standard applications. 
The examples presented here clearly show that a global constraint on sparsity, rather than quality 
constraints on the individual partition’s units, may beneht enormously the quality of the signal 
approximation. This is true even if the actual approximation of each unit is performed individually. 
Because considerations are restricted to producing high quality approximations, there is no risk that 
blocking artifacts could appear due to the signal division into non-overlapping pieces. 

Notational convention 

Throughout the paper M, C and N indicate the sets of real, complex, and natural numbers, respec¬ 
tively. Boldface letters are used to indicate Euclidean vectors, whilst standard mathematical fonts 
indicate components, e.g., f G iV G N is a vector of components f{i), i = 1,..., N. A partition 
of a signal f G is represented as a set of disjoint pieces i{q} G q = 1,... ,Q, which for sim¬ 
plicity are assumed to be all of the same size and such that QW = N. The signal is reconstructed 
from the partition through the operation f = where J represents the operator which con¬ 

catenates the partition. The operation is dehned as follows: given f{l} G and f{2} G 
the vector f = f{l} J f{2} is a vector in having components f{i) = for i = 1,..., W, 
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and f{i) = f{2}{i — Ny,) fori = Ny, + 1,..., 2N\,. Thus f = is a vector in having 

components f{i) = f{q}{i — {q — l)-^b), i = {q — l)-^b + 1, • • •, qNh, q = I,... ,Q. Consequently 
(f, f) = ||f|p = Ylq=i l|f{5'}|l^ ) were (•, •) indicates the Euclidean inner product and || ■ || the induced 
2-norm. 

Paper contributions 

Given a signal partition f{g} G q = 1,... ,Q and a dictionary V = |d„ G ; ||d„|| = to 

approximate the elements f{q} in the partition, the following outcome has been recently reported [15]: 
A very signihcant gain in the sparsity of a signal approximation may be effectively obtained by a 
greedy pursuit strategy, if the approximation of each piece f{g} (called ‘g-block’) is accomplished in a 
hierarchized manner. Suppose that the A;(g)-term approximation of f{5'} is the atomic decomposition: 

k(q) 

= X]c{9}nd£{q}„, g = l,...,Q, 

n=l 

with the atoms d(^{q}„, n = 1,..., A;(g) selected from the dictionary "D, via a stepwise greedy pursuit 
strategy. Suppose also that the number of atoms to approximate the whole signal f is a fixed value 
K, i.e., K = ^ possibility to handle this constraint is to consider a hierarchized selection 

of the pieces f{q} to be approximated in each approximation step. Some remarkable results of this 
strategy, which has been termed Hierarchized Block Wise (HBW) greedy strategy, are illustrated 
in [15] by approximating images using the greedy algorithms MB and OMP. When these methods 
are applied in the proposed HBW fashion are called HBW-MP and HBW-OMP, respectively. 

While [15] focusses on highlighting the suitability of the HBW-OMP/MP method when the image 
approximation is carried out in the wavelet domain, this Communication extends the method as well 
as the range of applicability. The extended techniques are shown to produce hugely sparse high 
quality approximation of music signals. This is accomplished with trigonometric dictionaries, which 
are endowed with the additional advantage of enhancing the competitiveness of the approach in terms 
of computational complexity. 


4 



The extension of the idea outlined in [15] comprises: 


• A revised version of the approach, which optimizes the ranking of the pieces fjo'}, q = 1,... ,Q 
for their step wise approximation without signihcant extra computational cost. Additionally, 
the selection of atoms is also optimized by including the Optimized Orthogonal Matching 
Pursuit (OOMP) criterion [16]. 

• A HBW backward approach for downgrading the approximation when required. This is realized 
in a stepwise optimized manner, by removing terms in the approximations of the selected blocks. 

• An alternative to the HBW strategy which consists of two stages. The hrst stage involves the 
approximation of the pieces f{g}, g = 1,..., Q, up to a tolerance error. The second stage rehnes 
the previous approximation by making possible the downgrading of the atomic decomposition of 
some blocks, called ‘donors’, and the upgrading of the atomic decompositions of another blocks, 
called ‘receivers’. Since this process is inspired in the Swapping-based-Rehnement of greedy 
strategies introduced in [17], we refer to it as HBW-Swapping-based-Rehnement (HBW-SbR). 

Further contributions of the paper are 

• The illustration of the huge gain in the sparsity of the representation of music signals obtainable 
by the proposed approach when using trigonometric dictionaries. 

• The hnding that a mixed trigonometric dictionary having both, discrete cosine and sine com¬ 
ponents, renders the best sparsity performance for melodic music signals, in comparison with 
the single discrete cosine (or sine) dictionaries with the same redundancy. This result contrasts 
with the fact that, the broadly used single discrete cosine basis is the basis yielding best spar¬ 
sity results for music signals, in comparison to other trigonometric basis of the same nature. 
Note: The discrete cosine/sine components of the mixed dictionary are not just the real and 
imaginary parts of the complex exponentials in a Fourier series. The involved phase factors 
make the whole difference in achieving high sparsity. 
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• The provision of implementation details in the form of algorithms for the proposed strategies 
and the provision of dedicated algorithms to operate with trigonometric dictionaries via the 
Fast Fourier Transform (FFT). 

A library of MATLAB functions for implementing the proposed methods, and running the numerical 
examples in this paper, are available for downloading on [18]. 

2 Hierarchized Blockwise OMP, revised 

As already stated, a signal f G will be considered to be the composition of Q identical and 
disjoint blocks f = where i{q} G assuming that N]^Q = N. Each piece i{q} G 

is approximated using a dictionary T) = {d„ G ; ||dn|| = with M > N^, by an atomic 

decomposition of the form; 

k{q) 

= 5^c{?}nd£{q}„, g = l...,g. (1) 

n=l 

where the atoms n = 1,..., k{q) are selected from the dictionary to approximate the block 

q. Thus, each set of selected indices F{g} = {^{q}n}^n=i is a subset of the set of labels which 

identify the atoms in V. In general F{g} ^ F{p} and k{q) ^ k{p) ii p ^ q. 

The HBW-OMP approach outlined in [15] selects atoms in (1) as indicated by the OMP approach, 
i.e.: On setting initially k{q) = 0 and r°{g} = f{g}, the algorithm picks the atoms in the atomic 
decomposition of block q by selecting one by one the indices which satisfy: 

k(q) 

^U}k(g)+i = argmax |(d„,r'^(‘^i{g})| , with r^^‘^'^{q} = i{q} - ^c{g}„d^{g} if k{q) > 0. (2) 

' ' n=l 

The coefficients c{g}^, n = 1,..., k{q) in (2) are such that the norm is minimized. This 

is ensured by requesting that = f{g} — Py-j i{q}, where Py? is the orthogonal projection 

k{q) k{q) 

operator onto = span This condition also guarantees that the set {d.£{q}n}n=i ^s 

linearly independent. 
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An implementation as in [16] provides us with two representations of One of the represen¬ 
tations is achieved by orthonormalization of the set The other by biorthogonalization 

of the same set. We implement the orthogonalization by Gram Schmidt method, including a re- 
orthogonalization step, as follows: The orthonormal set is inductively constructed from 

w{g}^ = w{g}^ = through the process: 

w{g}„ = d^{g}„-^w{g}.(w{g}.,d^{g}„>. (3) 

For numerical accuracy at least one re-orthogonalization step is usually needed. This implies to 
recalculate the above vectors as 

n—1 

^{Q}n ^ (w{g}i, w{g}„), (4) 

i=l 

with the corresponding normalization giving rise to the orthonormal set Notice that, 

whilst this set can be used to calculate the orthogonal projection of f{g} onto as 

k{q) 

(w{gL, f{g}), (5) 

n=l 

this superposition is not the atomic decomposition in terms of the selected atoms. In order to 
produce such a decomposition we use the other representation for an orthogonal projector given 
in [16]. Namely, 

k(q) 

(b{g}n, f{g}) • (6) 

n=l 

For a hxed q the vectors b{g}^, n = l,...,k{q) in (6) are biorthogonal to the selected atoms, 
i-e- (d<;{g}„, = 0 if n ^ m and 1 if n = m, and span the identical subspace i.e., = 

span{b{g}^}^^| = span{d£{g}„}^^j|. In order to fulhll the last condition all the vectors need to be 
updated when a new atom is introduced in the spanning set. Starting from b{g}^ = d^jg}^ both the 


calculation and upgrading of the biorthogonal set is attained recursively as follows 

as in (3) 

b{g}n ^ b{g}^ - b{g}^^ , n = 1,..., {k{q) - 1). 
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These vectors produce the required atomic decomposition for the block q. The coefficients in (1) are 
the inner products in (6), i.e. c{q}n = (b{g}„, , n = 1,..., k{q). 

The HBW variant of a greedy strategy establishes an order for upgrading the atomic decompo¬ 
sition of the blocks in the partition. Instead of completing the approximation of each block at once, 
the atomic decomposition to be upgraded at each iteration corresponds to a block, say q*, selected 
according to a greedy criterion. In other words: the HBW version of a pursuit method for approxi¬ 
mating a signal by partitioning involves two selection instances: a) the selection of dictionary’s atoms 
for approximating each of the blocks in the partition and b)the selection, at each iteration step, of 
the block where the npgrading of the approximation is to be realized. In [15] the blocks are selected 
throngh the process 


Q 


* 


arg max 

q=l,...,Q 




( 8 ) 


We notice, however, that condition (8) for the block’s selection can be optimized, in a stepwise sense, 
withont signihcant compntational cost. The next proposition revise the condition (8) considered 
in [15]. 


Proposition 1. Let i{q}be the index arising, for each value of q, from the maximization 
process (2). In order to minimize the square norm of the total residual at iteration k + 1 the 

atomic decomposition to he upgraded should correspond to the block q* such that 


g* = argmaxx(g), where xid) = (f{g} 




w{q} 


k(q)+l 


(9) 


with and as given in (3). 


Proof. Since at iteration k + 1 the atomic decomposition of only one block is upgraded by one atom, 
the total residne at iteration /c -|- 1 is constructed as 


r^+i = 

p¥=i 



Then, 


Q 

p=i 

p¥^<i 

Moreover, since = f{g} — using the orthonormal vectors (3) to calculate 

we have: 

which is minimum for the g*-value corresponding to the maximum value of f{g}^ . 

Since and Py.^^^ is hermitian (9) follows. □ 

Notice that, since the vectors ^{(l}k{q)+i will be used for subsequent approximations, the opti¬ 
mized ranking of blocks (9) does not require signihcant extra computational effort. Only Q — 1 of 
these vectors will not have been otherwise used when the algorithm stops. Apart from that, the 
complexity of selecting the blocks remains being that of hnding the maximum element of an array 
of length Q, i.e. 0{Q). 

Before presenting the algorithm details we would like to consider also the optimization of the 
selection of atoms (stage a) of a HBW pursuit strategy). Within our implementation the optimized 
selection of atoms is readily achievable by the Optimized Orthogonal Matching Pursuit (OOMP) 
approach [16], which selects the index i{q}k(q)+i through the maximization process: 

\(d ^ 

^{Q}kiq)+i = arg max — i , with s{q}n = ^ | (d„, w{g}J p. (10) 

n=l,...,M (l-s{g}n)2 

The OOMP selection criterion (10) minimizes, in a stepwise sense, the norm of the local error when 
approximating a single block [16]. The extra computational cost, in comparison to the OMP selection 
step (c.f. (2)) is the calculation of the sum s{q}n in the denominator of (10), for every atom in the 
dictionary, n = 1,..., M. As becomes clear in Algorithm 3, by saving the sum of previous iterations, 
at iteration k{q) + 1 only the inner products /d„, , n = 1,..., M need to be calculated. 
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2.1 Algorithmic Implementation 


The implementation details of the HBW-OOMP method is given in Algorithm 1, which is realized 
through Algorithms 2 and 3. For a hxed number K the algorithm iterates until the condition 
Ylq=i = K is met. 

2.2 Trigonometric dictionaries and sparse representation of music sig¬ 
nals 


The viability of a sparse representation for a given signal depends in large part on the dictionary 
choice. A possibility to secure a suitable dictionary for approximating a signal by partitioning 
is through learning techniques [19-22]. Another possibility is to build the dictionary by merging 
basis, or dictionaries, containing atoms of different nature. For instance, a Cosine-Dirac dictionary 
enlarged by the incorporation of B-spline dictionaries has been shown to be suitable for producing 
highly sparse high quality approximation of images [23]. Such dictionaries are very redundant. In 
the case of melodic music, however, stunning sparsity levels may be achieved with relatively much 
less redundant dictionaries. As will be illustrated here by numerical examples, it is the combination 
of a Redundant Discrete Cosine Dictionary (RDCD) and a Redundant Discrete Sine Dictionary 
(RDSD) D®, dehned below, which yields highly sparse representation of music signals when processed 
with the proposed techniques. 


i ^ cos V = 1 m,\m 


= \ ^ sin i = 1 


where w^{n) and w^{n), n = 1,... ,M are normalization factors as given by 


{ In. . 


if n = 1, 


iVb , M W -r / . 

2 + 2(1-cos(H^)) ^ 


y/Nh if n = 1, 

2 2(l-cos(^)) 


if n 7^ 1. 
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Algorithm 1 HBW-OOMP Procedure 

Input: Signal partition ^{q} G q = I,... ,Q. Dictionary V = {d„ G ; ||d„|| = 
Number K of total atoms to approximate the whole signal. 

Output: Sets r{q'}, q = containing the indices of the selected atom for each block. 

Orthonormal and biorthogonal sets and {h{q}^}'^^l\, for each block (c.f. (4) and (7) 

respectively). Coefficients of the atomic decompositions for each block q = 1,... ,Q. 

Approximation of the partition, P{g}, q = 1,... ,Q. Approximation of the whole signal P. 
{Initialization} 

j = 0 

for q = 1 : Q do 

k{q) = 1; r{g} = f{g}; s{q}n = 0, n = 1,..., M 
(Select the potential first atom for each block, as below} 
i{q} = argmax|(d„,f{g})| 

x{q) = |(d£{q}, fj^})! (Store the maximum for each block} 

Set w{g}^ = de{gy, w{g}^ = w{q}y b{g}^ = w{g}p P{g} = 0 

end for 

while j < K do 

{Select the block to be approximated (c.f.(9)) and store the index of the atom in the approxi¬ 
mation of that block, as below} 
g* = argmaxx(g) 

q=l,...,Q 

P{g*} ^ r{g*}U£{g*} 

(Update the residue corresponding to block q*, as below} 
r{g*} ^ r(g*} - w(g*}fc(g*) (w{g*}fc(g*), f{g*}) 
if k{q*) > 1 then 

(Upgrade the biorthogonal set {d{q*}n}n=i fo include the atom d^jg*} (c.f. (7))} 
end if 

(Apply Algorithm 2 to select the index i{q*} for a new potential atom for the block q*} 
Increase k{q*) ■(— k{q*) + 1 

(Compute w{g*}fc(g*) and w{g*}fe(g*) from d^^^*}, (c.f. (3) and (4)} 

(Update the objective functional y for the block q*} 

X{q*) = |(w{g*}fc(g*),f{g*})| 

Increase j j + 1 

end while 

(Calculation of coefficients and approximation} 
for q = 1 ■. Q do 

p{g} = f{g} - r{g}; c{q}n = (b{g}„, f{g}) , n = 1,..., k{q) 

end for 

(Concatenation of the partition to produce the whole signal approximation} 

f‘ = jjLif*{9} 
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Algorithm 2 Atom Selection Procedure 


Input: Residue r{q*}. Dictionary V = {d„ G ; ||d„|| = Auxiliary sequence 

{s{q*}n\n=i vector w{g*}fc(q*) for the upgrading. Set of already selected atoms Pjg*}. 
Output: Index £{g*} corresponding to the selected atom 

{Apply the Auxiliary Procedure, Algorithm 3, to upgrade the sequence {s{q*}n}n=i with respect 
to vector w{q*}kiq*)} 


(Select the index, as below} 

nr *1 l(dn, 

^{q } = argmax f 

n=l,...,M (1 - S{g*}n)2 
niV{q*} 



For M = Nt, each of the above dictionaries is an orthonormal basis, the Orthogonal Discrete Cosine 
Basis (ODCB) and the Orthogonal Discrete Sine Basis (ODSB), henceforth to be denoted and 
B^ respectively. The joint dictionary is an orthonormal basis for M = the Orthogonal Discrete 
Cosine-Sine Basis (ODCSB) to be indicated as For M > Ab, is a RDCD and P® a RDSD. 
If M > ^ UV^ becomes a Redundant Discrete Cosine-Sine Dictionary (RDCSD). 

For the sake of discussing a fast calculation of inner products with trigonometric atoms, given a 
vector y G M^, let’s dehne 

M 

.F(y,n,M) = n = l,...,M. (11) 

i=i 

When M = Ab (11) is the Discrete Fourier Transform of vector y G which can be evaluated 
using FFT. If M > Ab we can still calculate (11) via FFT by padding with zeros the vector y. Thus, 
(11) can be used to calculate inner products with the atoms in dictionaries D'’ and "D®. Indeed, 

^ 7r(2j_l^n—^ (^e“*^^J^(y,n,2M)j , n = 1,.. .,M. (12) 


12 








and 


^ 7r(2j_l^n—1) ^ j'g *" 2 VV(y,n, 2 M)j , n = 2,...,M + l, (13) 

1=1 

where Re(; 2 ) indicates the real part of and Ini( 2 ;) its imaginary part. 

The computation of inner products with trigonometric dictionaries via FFT is outlined in Algo¬ 
rithm 4. For dictionaries and that procedure reduces the complexity for calculating the inner 
products from O(iVbM), per block, to 0(2Mlog2 2M), per block. For dictionary the reduction 
is larger: 0(Mlog2M), because both the real and imaginary parts of the FFT are used. When 
Algorithm 1 is applied with trigonometric dictionaries the pieces implemented by Algorithms 2 and 
3 should be replaced by Algorithms 6 and 7, given in Appendix A. In addition to speeding the calcu¬ 
lations via FFT, Algorithms 6 and 7 avoid storing the dictionary. In the case of the dictionary it 
is assumed that both the Cosine and Sine components are of the same size and the hrst ^ elements 
of the dictionary are Cosine atoms. 

Algorithm 4 Computation of inner products with a trigonometric dictionary via FFT. 

IPTrgFFT procedure: p =IPTrgFFT(r, M, Case) 

Input: r G M, number of elements in the dictionary, and Case (‘Cos’ or ‘Sin’). 

Output: Vector p G with the inner products between r and ‘Cos’ or ‘Sin’ dictionaries. 

{Computation of auxiliary vector t G to compute p.} 

t = FFT(r, M) (FFT with (M — W) zero padding}. 

Case ‘Cos’ 

pW = TRTJ Re(e*^^t(n)), n = 1,..., M{(c.f. (12))} 

Case ‘Sin’ 

pW = Im(e*^^ t{n)), n = l,...,M {(c.f.(13))} 

Numerical Example I 

The signal approximated in this example is a piece of piano melody shown in Fig 1. It consists of 
N = 960512 samples (20 secs) divided into Q = 938 blocks with W = 1024 samples each. As sparsity 
measure for a signal representation we consider the Sparsity Ratio (SR) dehned as SR = where 
K is the total number of coefficients in the signal representation. As a measure of approximation 
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Figure 1: Piano melody. Credit: Julius O. Smith, Center for Computer Research in Music and Acoustics 
(CCRMA), Stanford University. 


quality we use the standard Signal to Noise Ratio (SNR), 

Ehdl/{9}(i)P 

q=l 

Afb,Q 
2 = 1 
q=l 

This numerical example aims at illustrating the following outcomes in relation to the signal in hand. 

1) The approximation power of all the orthogonal basis and dictionaries dehned above remarkably 
improve if, instead of applying OMP/OOMP independently to approximate each block f{g} 
up to the same SNR, the HBW-OMP/OOMP approach is applied to match the sparsity of the 
whole signal. 

2) For approximating the signal with the OMP/OOMP greedy strategies, the redundant dictio¬ 
naries and perform signihcantly better than any of the orthogonal basis and 

^cs 

In oder to demonstrate 1) and 2) each of the blocks i{q} is approximated independently with 
OMP/OOMP, up to SNR=25 dB. Redundant dictionaries, with redundancy 2 and 4, are simply 
creating by setting M = 2iVb and M = 4iVb in the dehnitions of and with iVb = 1024. They 
will be denoted as and respectively. Approximations of the same quality are per¬ 

formed with each of the ortho normal basis B^, B^ and 13“, and the mixed dictionaries = B^^UB^ 



SNR = 10 logi 


If - R 


= 10 logi 
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Diet. 

OMP 

HBW-OMP 

OOMP 

HBW-OOMP 


SR 

SNR 

SR 

SNR 

SR 

SNR 

SR 

SNR 


14.38 

25.0 

14.38 

35.19 

14.38 

25.0 

14.38 

35.19 


17.75 

25.0 

17.75 

34.91 

18.47 

25.0 

18.47 

35.08 

pc4 

19.39 

25.0 

19.39 

34.64 

19.80 

25.0 

19.80 

34.95 


7.65 

25.0 

7.65 

28.31 

7.65 

25.0 

7.65 

28.31 

ps2 

12.13 

25.0 

12.13 

29.39 

12.63 

25.0 

12.63 

29.39 

ps4 

13.17 

25.0 

13.17 

29.34 

13.42 

25.0 

13.42 

29.33 


10.77 

25.0 

10.77 

30.09 

10.77 

25.0 

10.77 

30.09 

•pcs2 

20.45 

25.0 

20.45 

35.81 

22.08 

25.0 

22.08 

35.56 


23.83 

25.0 

23.83 

35.56 

26.18 

25.0 

26.18 

36.37 


Table 1: Comparison of the approximation quality (SNR values) and sparsity (SR values) produced with 
trigonometric basis and redundant trigonometric dictionaries, and 2?“^. 

The second column shows the SR resulting when applying the block independent OMP approach to achieve 
a SNR=25.0dB. The fifth column demonstrates the significant gain in SNR rendered by the HBW-OMP 
strategy for the same sparsity. Further improvements of the same nature are shown in the last four columns 
corresponding to the approaches OOMP and HBW-OOMP. 


and U The results are presented in Table 1. As can be seen, the SR produced by 

the redundant dictionaries and is substantially larger than that corresponding to 

any of the orthogonal basis. Moreover, in all the cases the HBW-OMP/OOMP strategies improve 
notoriously (up to 11 dB) upon the OMP/OOMP approaches applied independently to produce a 
uniform SNR in the approximation of each block. 

It is noticed that the highest sparsity is attained by the dictionary which, as already discussed, 
involves the most effective implementation via FFT. It is also seen that HBW-OOMP over performs 
HBW-OMP in terms of sparsity (the running times are similar). This will appear even more clearly 
in Table 2, where the approximation of both methods are downgraded to produce a SNR of 25dB. 
The next section discusses the HBW backward strategy which allows for removing atoms from the 
signal approximation, in a stepwise optimized fashion. 

2.3 Hierarchized Blockwise Backwards OOMP 

We extend here the Backward Optimized Orthogonal Matching Pursuit (BOOMP) strategy [24] to 
select also the blocks from which the atoms are to be removed for downgrading the approximation 


15 




of a whole signal. The BOOMP strategy is stepwise optimal because it minimizes, at each step, 
the norm of the error resulting by downgrading, by one atom, an atomic decomposition. This 
result readily follows from the recursive equations for modifying vectors to account 

for the elimination of, say the j-th atom, from the set For each g, the reduced set of 

vectors spanning the reduced subspace = span{d^{g}^}^^| can be quickly obtained 

through the adaptive backward equations [24,25] 


h{q}n n = 1,..., j - 1, j + 1,..., k{q). (14) 

Consequently, the coefficients of the atomic decomposition corresponding to the block q, from which 
the atom is taken away, are modihed as 


c{9}.. <- - c{q}j - l.i + 1. ■".*:(«)■ 


(15) 


For a hxed value of q, the BOOMP criterion removes the coefficient c{q}jo such that [24] 


fU} 


Hq} 


argmm 


ji 


( 16 ) 


The equivalent criterion applies to the selection of the g*-block for the downgrading of the whole 
approximation. The proof of the next proposition parallels the proof of (16) given in [24]. 


Proposition 2. Assume that the approximation of a signal is given as P = where P{g} = 

Py?^ Let j^{q}, q = I,... ,Q be the indices which satisfies (16) for each q. The block q^, from 

where the atom is to be removed, in order to leave an approximation such that ||P —fl'H 

takes its minimum value, satisfies the condition: 


\(^{q}jA<i]\ 

q=l,...,Q l|b{g}jO{g}| 


q = argmm 


(17) 


Proof. Since at each step only one atom is removed from a particular block, say the g-th one, it holds 
that ||P — fj^ll = ||P{g} — The identical steps as in [24] lead to the expression 


P{q} - ifiq} 


d<( 




(b{gV|g},f{g}) 

l|b{g}go{g}IP 


(18) 
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2 


Thus IIPM-Wf = 


HQ}jO{q},i{q}) 


. Since (b{g}^.o{g}., f{(?}y = c{q}jo{qj it follows that 


I fa _ fa II jg jxiinimized by removing the atom corresponding to the block g* satisfying (17). □ 


The HBW-BOOMP algorithm for downgrading a given atomic decomposition of a signal is pre¬ 
sented in Algorithm 5. 

Numerical Example II 


In order to illustrate the backward approach, we change the information in the last four columns 
of Table 1 to have the sparsity of all the approximations corresponding to SNR=25dB. For this we 
downgrade the HBW-OMP/OOMP approximations to degrade the previous quality. As seen in Table 


Diet. 

OMP 

HBW-BOOMP 

OOMP 

HBW-BOOMP 


SR 

SNR 

SR 

SNR 

SR 

SNR 

SR 

SNR 

QC 

14.38 

25.0 

25.56 

25.0 

14.38 

25.0 

25.56 

25.0 

pc2 

17.60 

25.0 

34.05 

25.0 

18.47 

25.0 

36.34 

25.0 

pc4 

19.25 

25.0 

37.01 

25.0 

19.80 

25.0 

39.06 

25.0 

QS 

7.65 

25.0 

13.67 

25.0 

7.65 

25.0 

13.67 

25.0 

ps2 

12.03 

25.0 

25.74 

25.0 

12.63 

25.0 

27.56 

25.0 

ps4 

13.11 

25.0 

28.18 

25.0 

13.42 

25.0 

29.59 

25.0 

QCS 

10.77 

25.0 

19.94 

25.0 

10.77 

25.0 

19.94 

25.0 

-pcs2 

20.21 

25.0 

42.43 

25.0 

22.08 

25.0 

48.48 

25.0 

-pcs4 

23.53 

25.0 

50.19 

25.0 

26.18 

25.0 

59.68 

25.0 


Table 2: Comparison of sparsity (SR values) for the same SNR. The HBW-BOOMP results are obtained 
by degrading the HBW-OMP/OOMP approximation in Table 1 to SNR=25dB. 


2, the result is that the sparsity increases drastically. The SR improves up to 128% (for the 
dictionary) in comparison with the standard application of the same approaches without ranking the 
blocks. Notice also that the best SR result (for dictionary is 133% higher than the best SR 

result for an orthogonal basis (the basis). 

3 Refinement by Swaps 

The suitability of the HBW strategy for approximating the class of signals we are considering is a 
notable fact. However, the possibility of approximating every element of the partition completely in- 
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Algorithm 5 HBW-BOOMP 


Input: Biorthogonal sets q = Coefficients in the approximation of each 

block {c{q}n}t=\j Q = Sets of indices in the decomposition of each block r{g}, q = 

1,... ,Q. Number of total atoms for the downgraded approximation, K^. 

Output: Downgraded biorthogonal sets. Coefficients of the downgraded atomic decompositions 
for each block. Indices of the remaining atoms. 

{Initialization} 

A = 0; 

for q = 1 -. Q do 

k{q) = |r{g}| (Number of elements in r{g}} 

K = K + k{q) 

(Select the index of the potential hrst atom to be removed in each block, as below} 

■or T . \c{(l}n\ 

J {q} = argmm 

n=l,...Mq) l|b{g}nlr 


Xiiq) = 

end for 


\^{q}j^{q} I 

I 


(Store the minimum} 


i = K 

while i > Ad do 

(Select the block to be downgraded, as below} 

= argminxi(g) 

q=l,...,Q 

(Apply backward biorthogonalization (c.f. (14)) to downgrade the biorthogonal set correspond¬ 
ing to block q^, with respect to b{g^}jo{go} and shift indices} 

(Update the coefficients corresponding to block q^ (c.f. 15) and shift indices} 

Downgrade the set r{g^} r{g*}/t'{g^}jo{go} and shift indices 

Decrease k{q) A;(g) — 1 

(Select the new potential atom to be removed from block as below} 

■or oi • \c{q^}n\ 

3 [q }= argmm ^ 

n=l,...,k{q) ||b(g }n|| 

(Update the objective functional x; for the block g*, as below} 

’ WHq^Uq^r 

Decrease i ^ i — 1 


end while 





dependent of the others has the convenient feature of leaving room for straightforward parallelization. 
This is particularly important for implementations on Graphics Processing Units (GPU), for instance. 
Accordingly, an alternative worth exploring is to maintain the block independent approximation but 
allowing for possible eventual HBW rehnements. 

Assuming that the approximation (1) of each block in a signal partition is known, the goal now 
is to improve the signal approximation maintaining unaltered the total number of atoms K. The 
proposed rehnement consists of movement of atoms controlled by the following instructions: 

i) Use criterion (17) to remove one atom from the block g*. Let us call this block a ‘donor’ block 
and indicate it with the index gd- 

ii) Use criterion (9) to incorporate one atom in the approximation of the block g*. Let us call this 
block a ‘receiver’ and indicate it with the index q,.. 

iii) Denote by (5d the square norm of the error introduced at step i)by downgrading the approx¬ 
imation of the donor block gd. Denote by the square norm of the gain in improving the 
approximation of the receiver g^. Proceed according to the following rule: 

If 6t. > 5d accept the change and repeat steps i) and ii). Otherwise stop. 

The above procedure is termed HBW-SbR-OMP/OOMP according to whether the selection of the 
atom at stage ii) is realized through the OMP (c.f. (2)) or OOMP (c.f. (10)) criterion. 

Notice that the HBW-SbR-OMP/OOMP methods, implemented as proposed above, possess the 
desirable feature of acting only if the approximation can be improved by the proposed swap. The 
possible situation for which g^ = gd corresponds to a swap in the original SbR-OMP/OOMP approach 
[17], which may take place here as a particular case. 

Implementation Remark: For successive implementation of the steps i) and ii) one needs: to 
downgrade the vectors b{g*}„, n = 1,... k{q) at step i) and to upgrade these vectors at step ii). In 
order to allow for step ii) the orthogonal set has to be downgraded as well. This can 
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be realized effectively by the plane rotation approach [26,27]. 

Numerical Example III 


We illustrate the HBW-SbR-OMP/OOMP approaches by applying them to the signal in Fig. 2, 
which is a flute exercise consisting of iV = 96256 samples divided into Q = 94 blocks with W = 1024 
samples each. 

Table 3 and 4 show the improvement in quality obtained by applying the HBW-SbR-OMP/OOMP 



01 23456789 10 

X 10“ 


Figure 2: Flute exercise: N = 96256 samples. 


approaches on outputs of the OMP/OOMP approximations. Notice that in most cases the SNR 
results are practically equivalent to those obtained by the HBW-OMP/OOMP approaches. 


Diet. 

SR 

OMP 

Swaps 

HBW-SbR-OMP 

HBW-OMP 

QC 

10.13 

25.0 

1426 

26.80 

26.79 

pc2 

14.26 

25.0 

1013 

26.69 

26.67 

pc4 

15.58 

25.0 

945 

26.63 

26.61 

QS 

7.27 

25.0 

2612 

26.72 

26.73 

ps2 

12.10 

25.0 

1397 

27.24 

27.27 

ps4 

13.16 

25.0 

1435 

27.27 

27.25 

^cs 

8.08 

25.0 

1822 

26.65 

26.65 


29.00 

25.0 

636 

27.38 

27.37 

-pcs4 

35.17 

25.0 

518 

27.36 

27.36 


Table 3: Comparison of quality (SRN values) for a fixed sparsity: that corresponding to SRN=25 with 
the OMP approach. The forth column shows the number of swaps and the fifth columns the SRN achieved 
by those swaps through the HBW-SbR-OMP refinement to the outputs yielding the second column. For 
further comparison the last column shows the results corresponding to the HBW-OMP approach. 
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Diet. 

SR 

OOMP 

Swaps 

HBW-SbR-OOMP 

HBW-OOMP 

QC 

10.13 

25.0 

1426 

26.80 

26.79 


14.83 

25.0 

865 

26.66 

26.68 

pc4 

15.80 

25.0 

929 

26.60 

26.61 


7.27 

25.0 

2612 

26.72 

26.73 

ps2 

12.58 

25.0 

1531 

27.14 

27.27 

ps4 

13.34 

25.0 

1445 

27.27 

27.26 


8.08 

25.0 

1822 

26.65 

26.65 

-pcs2 

33.47 

25.0 

556 

27.64 

27.63 


42.05 

25.0 

228 

27.10 

27.69 


Table 4: Same description as in Table 3 but for the OOMP, HBW-SbR-OOMP, and HBW-OOMP ap¬ 
proaches. 


4 Processing of large signals 

As already discussed, the additional computational complexity introduced by the HBW raking of 
Q blocks is only 0((5) (over the complexity of a standard pursuit algorithm implemented without 
ranking the blocks). Storage does become more demanding though. The need for saving vectors (7) 
and (3), for instance, implies to have to store 2K vectors of size Af,. This requirement restricts the 
amount of blocks to be processed with a small system such as a standard laptop. Hence, we outline 
here a processing scheme which makes it possible the application of HBW techniques on very large 
signals. For this the whole signal needs to be divided into large segments of convenient size, say 
Ng. More specihcally, a number of, say P blocks, of size W are grouped together to form S larger 
segments g{s} = s = 1,..., A, where it is assumed that = Ng and SNg = N\^Q. 

Each segment g{<s} G is approximated with a HBW strategy as an entire individual signal. The 
approximation of the original signal is assembled at the very end of the process by the operation 

^ S 

P = In other words, the processing of large signals is realized by chopping the signal 

into segments and processing each segment independently of the others. Nevertheless, the question 
as to how to set the sparsity constraint for each segment needs further consideration. One could, 
of course, require the same sparsity in every segment, unless information advising otherwise were 
available. While uniform sparsity guarantees the same sparsity on the whole signal, in particular 
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cases, where the nature of the signal changes over its range of dehnition, it would not render the 
best approximation. In order to avoid encountering this type of situation we introduce a preliminary 
step: the randomization of all the small blocks in the signal partition. The implementation is carried 
out as follows: 

i) Given a signal f split it into Q blocks of size N\^. Apply an invertible random permutation 11 
to scramble the block’s location (the bottom graph in Figure 3 provides a visual illustration of 
how the signal in the top graph looks after this step). Let’s denote the re-arranged signal by 

f = AiG}. 

^gp 

ii) Group every P of the blocks f{g}, q = 1,..., Q to have S segments g{s} = Jq=(s_i)p+if{5'}, s = 

1,...,S = 2 

hi) Approximate each of the segments g{<s} independently of the others, as it each segment were 

'^SP 

an entire signal, i.e., obtain g’^js} = s = 1,..., S'. 

iv) Assemble the segments to have the approximation of the whole signal as f®' = 

v) Reverse the permutation of the block’s location to obtain, from P, the approximation of the 
original signal 



Figure 3: The top graph shows a piece (1433600 samples) of Piazzola music interpreted by Gidon Kremmer 
and orchestra. The bottom graph shows the resulting signal after the randomization of the blocks. 


22 




Numerical Example IV 

The signal to be approximated is a 32.5 secs (1433600 samples) piece of Piazzola music, shown 
in the top graph of Figure 3. It is partitioned into Q = 1400 blocks of 1024 points each. Af¬ 
ter randomization (bottom graph of the same figure) the blocks are grouped into 50 segments, 
each of which is independently approximated by the HBW-OMP/OOMP approach to produce a 
SR=11.53. The resulting SNR of the approximated signal is 29.12 dB with HBW-OMP and 30.15 dB 
with HBW-OOMP. Equivalent approximation quality is obtained with other realizations of the ran¬ 
dom process. For appreciation of this result the original signal was also approximated by HBW- 
OMP/OOMP, but without segmentation. In that case the quality corresponding to SR=11.53 is 
only 1.9% higher (SNR=29.66 dB with HBW-OMP and SNR=30.73 dB with HBW-OOMP). The 
conclusion is that, even if the largest the segments the less the distortion for the same sparsity, apply¬ 
ing HBW-OMP/OOMP on segments of a large signal may be still signihcantly more advantageous 
than the independent approximation of the blocks. Indeed, in this example the latter yields a SNR 
of 25 dB with OMP and 26.1 dB with OOMP. 

5 Conclusions 

Cooperative greedy pursuit strategies for approximating a signal partition subjected to a global 
constraint on sparsity have been considered. The cooperation between partition units was real¬ 
ized in several ways: i)By ranking the partition units for their sequential stepwise approximation 
(HBW-OMP/OOMP) ii)By ranking the partition units for stepwise downgrading of the whole ap¬ 
proximation (HBW-BOOMP) iii)By allowing for downgrading of some partition units (donors) and 
upgrading of another partitions units (receivers), to improve the approximation quality (HBW-SbR- 
OMP/OOMP). Comparisons of the OMP and OOMP criteria indicate that for the type of dictionaries 
and signals considered here, the OOMP criterion renders improvements of sparsity up to 20% in re¬ 
lation to the equivalent strategy involving OMP. 
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The HBW algorithms maintain the complexity comparable with the totally independent approxi¬ 
mation of each partition unit. The examples presented here, using trigonometric dictionaries for 
achieving high quality approximation of musics signals, provide a clear illustration of the benehts 
arising by the proposed strategies. Indeed, for a piano melody the increment in sparsity is 128% 
higher than the sparsity produced by approximating each partition unit up to same quality. Com¬ 
parison of trigonometric dictionaries results, against results obtained with trigonometric basis, give 
an improvement of equivalent order: 133%. The mixed Cosine-Sine dictionary has yielded the best 
sparsity results, not only for the test signals included in this paper but for a signihcant number of 
other signals, all in the category of melodic music. 

Because the proposed HBW implementation of pursuit strategies requires storage of order N\,K (W 
being the length of each block in the partition and K the number of total atoms to approximate the 
whole signals) a simple procedure for processing large signal was proposed; Previous to dividing the 
signal into independent segments, each of which to be approximated with the identical constraint of 
sparsity, the randomization of the partition units is advised. 
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Appendix A 


Algorithm 6 Atom Selection Procedure via FFT 


Input: Residue rjg*}. Dictionary Case (‘Cos’, ‘Sin’ or mixed ‘Cos-Sin’) and number of elements, 
M. Sequence {s{q*}n}n=i if Case is ‘Cos’ or ‘Sin’, and {s^{q*}n}n=ii {^^{Q*}n}n=i if Case is ‘Cos- 
Sin’. Vector to upgrade the sequences. Set of already selected atoms F{g*} (sets F^’lg*} 

and with the atoms in each component if Case is ‘Cos-Sin’). 

Output: Selected index i{q*} for a potential new atom for block {g*}. 

Cases single ‘Cos’ or ‘Sin’ 

{Call IPTrgFFT procedure. Algorithm 4, to calculate inner products, and apply Algorithm 7 to 
upgrade the sequence {s{q*}n}n=i with respect to vector w{g*}fc(g*)} 
p=IPTrgFFT(r{g*}, 2M, Case), 

(Select the index for a new potential atom for the block q*, as below.} 

i{q*} = argmax 

n=l,...,M (1 - S{q*}n)^ 
n^r{q*} 

Case ‘Cos-Sin’ 

(Same procedure as for the previous case but for both ‘Cos’ and ‘Sin’ cases}, 
p'’{g*}=IPTrgFFT(r{g*}, M ‘Cos’); p^{g*}=IPTrgFFT(r{g*}, M, ‘Sin’) 

(Apply Algorithm 7 to upgrade the sequences {s‘^{q*}n}n=i {^^{Q*}n}n=i selec the index} 

M 


Set M ^ 

i {q } = argmax^-| = argmax 


=1,...,M (1 - s®(g*}n)2 
n^rqg*} 


{1 - s^{q*}e^{q*})2 


a=l,...,M (1 - S'={g*}n)2 
(Evaluate the maximum value, as below} 

P \q s , , , 1) h \y j 

(1 (9 }<?c{g*})2 

/i{g*} = max(/i‘’(g*},/i'^(g*}) 
if /i = /i'’ theu 
= t{q*} 

else 


7(g*} = C{g*} + M 

eud if 
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Algorithm 7 Auxiliary Procedure via FFT 

Input: Vector ^{q*}k{q*)- Dictionary Case: ‘Cos’, ‘Sin’ or ‘Cos-Sin’, and number of elements 
M. Sequence {s{q^}n}n=ii Cases ‘Cos’ and ‘Sin’, or {s'^{q^}n}n=n {s^{q*}n}n=i^ C!ase 

‘Cos-Sin’. 

Output: Upgraded sequence, {s{q*}n}n=ij Cases ‘Cos’ and ‘Sin’, or {s^{q*}n}n=i 
for Case ‘Cos-Sin’. 

Case ‘Cos’ or ‘Sin’ 

{Call IPTrgFFT procedure. Algorithm 4, to calculate inner products with vector ^{q*}k{q*)} 
p{g*}=IPTrgFFT(w{g*}fc(q*), 2M, Case) 

(Upgrade the sequence, as below} 
for n = 1 : M do 

s{q*}n = s{q*}n+\p{q*}{n)f 
end for 
Case ‘Cos-Sin’ 

(Call IPTrgFFT procedure. Algorithm 4, to calculate the inner products between vector w{q*}k{q*) 
and each of the dictionary components: Cases ‘Cos’ and ‘Sin’} 
p^=IPTrgFFT(w{g*}fc(,*),M, ‘Cos’) 
p*=IPTrgFFT(w{g*}fc(,*),M, ‘Sin’) 

(Upgrade the sequences, as below} 
for n = 1 : ^ do 

sHq*}n = s^{q*}n + \pHq*}{n)f 
s^{q*}n = s^{q*}n + \f{q*}{n)f 
end for 
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